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(9) Summary of Research: 


My goal is to develop and implement efficient, accurate, and robust IMplicit-EXplicit Runge- 
Kutta (IMEX RK) methods [9] for overcoming geometry-induced stiffness with applications to 
computational electromagnetics (CEM), computational fluid dynamics (CFD) and computational 
aeroacoustics (CAA). IMEX algorithms solve the non-stiff portions of the domain using explicit 
methods, and isolate and solve the more expensive stiff portions using implicit methods. 

Current algorithms in CEM can only simulate purely harmonic (up to 10GHz plane wave) EM 
scattering by fighter aircraft, which are assumed to be pure metallic shells, and cannot handle the 
inclusion of coatings, penetration into and radiation out of the aircraft. Efficient IMEX RK 
methods could potentially increase current CEM capabilities by 1-2 orders of magnitude, 
allowing scientists and engineers to attack more challenging and realistic problems. 

This year, I completed my third year of research under the guidance of Professors David Gottlieb 
and Jan S. Hesthaven of Brown University and Dr. Mark H. Carpenter of NASA Langley 
Research Center. During the past 3 years, I implemented and tested explicit, implicit, and IMEX 
time-integration algorithms for solving linear as well as nonlinear equations in one and two 
dimensions on unstructured grids, such as burgers equation and the Euler equations, using 
Discontinuous Galerkin [3, 4, 5, 8] spectral element spatial discretizations. 

The Discontinuous Galerkin Spectral Element (DGSE) method builds upon the strengths and 
overcomes the weaknesses of the Discontinuous Galerkin Finite Element method [3, 4] and the 
classical spectral element method introduced by Patera [11, 12], and has a number of advantages 
over classical finite difference and finite volume methods. DGSE methods are especially well 
suited for IMEX algorithms, since they allow for clean and easy decoupling of the stiff from the 
nonstiff regions of the domain. Furthermore, they are highly parallelizable and accurate, provide 
for simple treatment of boundary conditions, handle complicated geometries well, and can easily 
handle adaptivity. Utilizing an unstructured grid [8] allows me to capture very fine details of 
complex domains, which is crucial for tackling real-world problems. 

My research shows that when geometry-induced stiffness is significant (greater or equal to 2 
orders of magnitude), IMEX algorithms [9] outperform traditional explicit time-stepping RK 
methods by about an order of magnitude or more in ID for systems having smooth solutions. I 
also found that Krylov subspace iterative methods, such as GMRES and BiCGStab, in 
conjunction with good preconditioners are essential in order to achieve this kind of speedup. 
Some of the results from my research have been presented at the Second MIT conference on 
Computational Fluid and Solid Mechanics and are published in [2], 

In 2D, IMEX methods outperform explicit methods even without preconditioning (when 
geometry-induced stiffness is greater than 1-2 orders of magnitude) for problems having smooth 
solutions. When shocks are present, the situation is not as clear since Newton’s method often 
fails to converge and many steps must be repeated with smaller time-steps. Preconditioning gives 
IMEX methods an even greater computational advantage in two dimensions. ILU 
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preconditioners worked very well in ID, while Jacobi, Gauss-Seidel and block-Jacobi 
preconditioners were not as effective. I will need to further investigate effective 2D and 3D 
preconditioning techniques such as multi grid-based preconditioners [10]. 

I am currently developing and testing time advancement algorithms for solving systems of 
nonlinear partial differential equations in multiple dimensions, such as Eulers equations and the 
compressible Navier-Stokes equations, using Discontinuous Galerkin [8] spectral element spatial 
discretizations. 


2D Nozzle Flow and Conclusion: 


The three-dimensional Navier-Stokes (NS) equations describe the behavior of many types of 
fluids and are given below: 
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Re = = Reynolds number, Pr = E = p(T+ 4(« 2 + r- + w 2 )) = total energy, 

u = dynamic viscosity, A = bulk viscosity, k = coefficient of thermal conductivity, 

stress tensor elements: r.,, = + Sy A s u = Kronecker delta, 

— ' > _ heat capacity (constant pressure 
1 ’ heat capacity (constant volume ) 


For our test case, we solve the two-dimensional Euler equations (in a nozzle), which can be 
derived from the NS equations by taking the limit of Reynolds number going to infinity (right- 
hand side in above equations becomes 0). 
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The axisymmetric converging-diverging nozzle has area given by: 

Area(x) = 1.75 - .75 * cos( (.2 * x- 1.0) * pi ), 0sxs5.0 

Area(x) = 1 .25 - .25 * cos( (.2 * x - 1.0) * pi ), 5.0 ^ x ^ 10.0 

We solve a steady-state (SS) problem where the SS solution has a shock at x = 7.56. The inflow 
Mach number = .240 and the outflow Mach number = .501 . The test case until the physical time t 
= 5, which is well before the time the shock develops. So for our test, the solution is relatively 
smooth and does not contain shocks, although it does contain some complex structures. 


Density 



x— Momentum 




Figure 1: IMEX, N-8, No Preconditioning, Newton Tot. - 1e-3, Bicgstab, Failed-O, nozzles 2, 
CPU-time- 11,1 15s (5,194 time steps) 


The implicit elements are located within a semicircular region of radius r = .25 centered at x = 
7.562 (center lies on bottom wall/centerline, see figure). The geometric stillness, which is 
defined as the ratio of average IMEX time step to average explicit time step is approximately 
2.27 (N = 8) for the unstructured triangular mesh used. The ratio of implicit to total elements is 
56/496, which is slightly greater than 10%. 

We use the method of lines to discretize the partial differential equations in space and time. First, 
we discretize space using a nodal Discontinuous Galerkin Spectral Element method based on [8]. 
Next, we integrate the system of ordinary differential equations (ODEs) in time. For the IMEX 
method, we use a modified Newton-Krylov method (Newton tolerance = IE-03) to linearize the 
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nonlinear system of ODEs, and to iteratively solve the linear systems (BiCGStab). Both time- 
integration methods are globally 4 th -order accurate. Preconditioners were not used. 

We impose boundary conditions on the characteristic variables at the inflow and outflow, and 
penalize the velocity against its mirror image at the top and bottom walls. 

We use a stability time-step controller for all runs: 


At = CFL/(A*N 2 *geometric factor), where A is the maximum wave speed. 


We can see that for N = 8 both methods take approximately the same CPU-time, with the IMEX 
being slightly faster. For N = 4, the explicit method is faster. 


Time- 

Integration 

Method 

Polynomial 
Order (N) 

Average 

Time-Step 

(At) 

Total Time- 
Steps 

Total CPU- 
time (sec.) 

Purely Explicit 

4 

1.39E-03 

3,601 

1,349 
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4.24E-04 

11,800 

11,313 

IMEX 

4 

3.34E-03 

1,499 

1,864 
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9.63E-04 

5,194 

11,115 


IMEX methods are not faster in this particular test, but will be significantly faster for the 
following situation: 


(1) Ratio of implicit to total elements < 5-10%. 

(2) Effective Preconditioners. 

(3) Geometric Stiffness > 2. 

(4) Navier Stokes equations, which have viscous terms. 

(5) Advanced Time-Step Controller (PDD). 


We believe that when some or all of the above are implemented, IMEX methods will become 
significantly more efficient than purely explicit time-integration methods. 
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